##############################################################
##############################################################
####### Replication code for "Do museums promote reconciliation? Evidence from a field experiment," Journal of Politics
####### This file includes code for all analyses and figures in the main text
####### Plots based on those used by Broockman and Kalla (2016)
##############################################################
##############################################################

require("sandwich")
require("plyr")
require("lmtest")
require(dplyr)
require(gridExtra)
require("RColorBrewer")
require(ggplot2)

##############################################################
###### Read in data
##############################################################
load(file = "all.Rdata")
##############################################################
###### Establish main functions
##############################################################

### ATE FUNCTIONS ##
# This estimates ATE when we have a pre-treatment measurement
est.ate<-function(dv, predv, df){
  predv <- f(predv)
  dv <- f(dv)
  summary(fit.1 <- lm(dv~df$treat + df$pre_ideology_1 +
                        + predv*df$date_diff + df$base_gender +df$age + df$v))
  vcv <- vcovHC(fit.1)
  n <- nobs(fit.1)
  result <- coeftest(fit.1, vcv)[2, 1:4] / itt.d
  result <- list("point.estimate" = result[1],"se"=result[2],"pvalue"=result[4], "obs" = n)
  return(result)
}
# This estimates ATE when we don't have a pre-treatment measurement
est.ate.np<-function(dv, df){
  dv <- f(dv)
  summary(fit.1 <- lm(dv~df$treat + df$pre_ideology_1 + df$base_gender +df$age+df$v))
  vcv <- vcovHC(fit.1)
  n <- nobs(fit.1)
  result <- coeftest(fit.1, vcv)[2, 1:4] / itt.d
  result <- list("point.estimate" = result[1],"se"=result[2],"pvalue"=result[4], "obs" = n)
  return(result)
}

# this function recodes NAs to the mean, per our PAP
f <- function(x){ 
  m <- mean(x, na.rm = TRUE) 
  x[is.na(x)] <- m 
  x 
}
## recode covariates to means
all$age <- f(all$age)
all$pre_ideology_1 <- f(all$pre_ideology_1)
all$base_gender <- f(all$base_gender)
all$date_diff <- f(all$date_diff)
# split dataset into left, right, and related to victim for heterogeneous analyses
left <- all[all$right == 0,]
right <- all[all$right == 1,]
itt.d <- all$itt.d
##############################################################
###### Figure 1: Political institutions treatment effects 
##############################################################

# Overall
dem <- est.ate(all$democracy, all$pre_democracy, all)
mil <- est.ate(all$military_gov, all$pre_military_gov, all)
gov_sat <- est.ate(all$inst_gov, all$pre_inst_gov, all)
mil_sat <- est.ate(all$inst_mil, all$pre_inst_mil, all)
pol_sat <- est.ate(all$inst_police, all$pre_inst_police, all)
gov_trust <- est.ate(all$conf_gov, all$pre_conf_gov, all)
mil_trust <- est.ate(all$conf_mil, all$pre_conf_mil, all)
pol_trust <- est.ate(all$conf_police, all$pre_conf_police, all)
church_trust <- est.ate(all$conf_church, all$pre_conf_church, all)

# right and left
dem.right <- est.ate(right$democracy, right$pre_democracy, right)
mil.right <- est.ate(right$military_gov, right$pre_military_gov, right)
gov_sat.right <- est.ate(right$inst_gov, right$pre_inst_gov, right)
mil_sat.right <- est.ate(right$inst_mil, right$pre_inst_mil, right)
pol_sat.right <- est.ate(right$inst_police, right$pre_inst_police, right)
gov_trust.right <- est.ate(right$conf_gov, right$pre_conf_gov, right)
mil_trust.right <- est.ate(right$conf_mil, right$pre_conf_mil, right)
pol_trust.right <- est.ate(right$conf_police, right$pre_conf_police, right)
church_trust.right <- est.ate(right$conf_church, right$pre_conf_church, right)

dem.left <- est.ate(left$democracy, left$pre_democracy, left)
mil.left <- est.ate(left$military_gov, left$pre_military_gov, left)
gov_sat.left <- est.ate(left$inst_gov, left$pre_inst_gov, left)
mil_sat.left <- est.ate(left$inst_mil, left$pre_inst_mil, left)
pol_sat.left <- est.ate(left$inst_police, left$pre_inst_police, left)
gov_trust.left <- est.ate(left$conf_gov, left$pre_conf_gov, left)
mil_trust.left <- est.ate(left$conf_mil, left$pre_conf_mil, left)
pol_trust.left <- est.ate(left$conf_police, left$pre_conf_police, left)
church_trust.left <- est.ate(left$conf_church, left$pre_conf_church, left)

# make figure 
results.df <- as.data.frame(rbind(dem[1], dem.left[1], dem.right[1],
                                        mil[1],  mil.left[1], mil.right[1],
                                        gov_sat[1], gov_sat.left[1], gov_sat.right[1],
                                        mil_sat[1],  mil_sat.left[1], mil_sat.right[1],
                                        pol_sat[1],  pol_sat.left[1], pol_sat.right[1],
                                        gov_trust[1],  gov_trust.left[1],gov_trust.right[1],
                                        mil_trust[1],  mil_trust.left[1], mil_trust.right[1],
                                        pol_trust[1],  pol_trust.left[1], pol_trust.right[1],
                                        church_trust[1],  church_trust.left[1], church_trust.right[1]
                                        ))
results.df$se <- c(dem[2], dem.left[2], dem.right[2],
                   mil[2],  mil.left[2], mil.right[2],
                   gov_sat[2], gov_sat.left[2], gov_sat.right[2],
                   mil_sat[2],  mil_sat.left[2], mil_sat.right[2],
                   pol_sat[2],  pol_sat.left[2], pol_sat.right[2],
                   gov_trust[2],  gov_trust.left[2],gov_trust.right[2],
                   mil_trust[2],  mil_trust.left[2], mil_trust.right[2],
                   pol_trust[2],  pol_trust.left[2], pol_trust.right[2],
                   church_trust[2],  church_trust.left[2], church_trust.right[2])
results.df$se <- unlist(results.df$se)
results.df$point.estimate <- unlist(results.df$point.estimate)
results.df$Variable <- NA
results.df$xpos <- NA
for (i in 1:3){results.df$Variable[i] <- "democracy"}
for (i in 4:6){results.df$Variable[i] <- "military govt"}
for (i in 7:9){results.df$Variable[i] <- "govt satisfaction"}
for (i in 10:12){results.df$Variable[i] <- "military satisfaction"}
for (i in 13:25){results.df$Variable[i] <- "police satisfaction"}
for (i in 16:18){results.df$Variable[i] <- "govt trust"}
for (i in 19:21){results.df$Variable[i] <- "military trust"}
for (i in 22:24){results.df$Variable[i] <- "police trust"}
for (i in 25:27){results.df$Variable[i] <- "church trust"}

results.df$varnum<- with(results.df, paste0(as.numeric(factor(Variable))))
results.df$varnum <- as.numeric(results.df$varnum)
results.df$sample <- rep(1:3, 9)
results.df$Population <- mapvalues(results.df$sample, c(1,2,3), c("All", "Left", "Right"), warn_missing = TRUE)
results.df$xpos<- as.numeric(paste(results.df$varnum, results.df$sample, sep = "."))
results.df$ci.hi <- unlist(results.df$point.estimate) + unlist(results.df$se) * 1.96
results.df$ci.low <- unlist(results.df$point.estimate) - unlist(results.df$se) * 1.96
results.df$se.hi <- unlist(results.df$point.estimate) + unlist(results.df$se)
results.df$se.low <- unlist(results.df$point.estimate) - unlist(results.df$se)
labels <- unique(results.df$Variable)
g <- ggplot(results.df,
            aes(x=xpos, y=point.estimate,
                group=Variable, color=Population)) +
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  geom_linerange(aes(ymin=ci.low, ymax=ci.hi)) +
  geom_linerange(aes(ymin=se.low, ymax=se.hi),lwd=1) +
  geom_point(color="black") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  ylab("Coefficient") +
  xlab("Variable") +
  ggtitle("Political Institutions")
g +  scale_color_grey() +
  ggplot2::annotate("text",
           label = c("church trust \n ctrl mean=1.15", "dem satisfaction \n ctrl mean=1.16", "govt satisfaction \n ctrl mean=0.80",
                     "govt trust \n ctrl mean=0.98", "mil govt (0-1) \n ctrl mean=0.28", "mil satisfaction \n ctrl mean=1.18",
                     "mil trust \n ctrl mean=1.17", "pol satisfaction \n ctrl mean=1.56", "pol trust \n ctrl mean=1.67"),
           x = c(1:9)+.2, y = -.5,
           colour = "black", size = 2.8) + theme_bw()
ggsave("polinst.pdf", width = 10, height = 6)

##############################################################
###### Figure 2: Transitional justice treatment effects 
##############################################################

advance <- est.ate.np(all$justice_advance, all)
justice_account <- est.ate.np(all$justice_account, all)
compensation <- est.ate(all$current_recomp, all$pre_current_recomp, all)
judicial <- est.ate.np(all$tj_judicial, all)
inst_apology <- est.ate.np(all$tj_apology, all)
apologize <- est.ate.np(all$policies_apologize, all)
compensate <- est.ate.np(all$policies_compensate, all)
pardoned <- est.ate.np(all$policies_pardon, all)

advance.right <- est.ate.np(right$justice_advance, right)
justice_account.right <- est.ate.np(right$justice_account, right)
compensation.right <- est.ate(right$current_recomp, right$pre_current_recomp, right)
judicial.right <- est.ate.np(right$tj_judicial, right)
inst_apology.right <- est.ate.np(right$tj_apology, right)
apologize.right <- est.ate.np(right$policies_apologize, right)
compensate.right <- est.ate.np(right$policies_compensate, right)
pardoned.right <- est.ate.np(right$policies_pardon, right)

advance.left <- est.ate.np(left$justice_advance, left)
justice_account.left <- est.ate.np(left$justice_account, left)
compensation.left <- est.ate(left$current_recomp, left$pre_current_recomp, left)
judicial.left <- est.ate.np(left$tj_judicial, left)
inst_apology.left <- est.ate.np(left$tj_apology, left)
apologize.left <- est.ate.np(left$policies_apologize, left)
compensate.left <- est.ate.np(left$policies_compensate, left)
pardoned.left <- est.ate.np(left$policies_pardon, left)

results.df.tj <- as.data.frame(rbind(advance[1], advance.left[1],advance.right[1],
                                           justice_account[1], justice_account.left[1], justice_account.right[1],
                                           compensation[1],  compensation.left[1], compensation.right[1],
                                           judicial[1],judicial.left[1], judicial.right[1],
                                           inst_apology[1],  inst_apology.left[1], inst_apology.right[1],
                                          apologize[1],  apologize.left[1], apologize.right[1],
                                           compensate[1],  compensate.left[1], compensate.right[1],
                                           pardoned[1], pardoned.left[1],pardoned.right[1]))

results.df.tj$se <- as.data.frame(rbind(advance[2], advance.left[2],advance.right[2],
                                     justice_account[2], justice_account.left[2], justice_account.right[2],
                                     compensation[2],  compensation.left[2], compensation.right[2],
                                     judicial[2],judicial.left[2], judicial.right[2],
                                     inst_apology[2],  inst_apology.left[2], inst_apology.right[2],
                                     apologize[2],  apologize.left[2], apologize.right[2],
                                     compensate[2],  compensate.left[2], compensate.right[2],
                                     pardoned[2], pardoned.left[2],pardoned.right[2]))
results.df.tj$se <- unlist(results.df.tj$se)
results.df.tj$point.estimate <- unlist(results.df.tj$point.estimate)
results.df.tj$Variable <- NA
results.df.tj$xpos <- NA
for (i in 1:3){results.df.tj$Variable[i] <- "advance"}
for (i in 4:6){results.df.tj$Variable[i] <- "accountable"}
for (i in 7:9){results.df.tj$Variable[i] <- "compensation"}
for (i in 10:12){results.df.tj$Variable[i] <- "punish"}
for (i in 13:15){results.df.tj$Variable[i] <- "public apology"}
for (i in 16:18){results.df.tj$Variable[i] <- "forced apology"}
for (i in 19:21){results.df.tj$Variable[i] <- "forced compensation"}
for (i in 22:34){results.df.tj$Variable[i] <- "pardoned"}

# X position of different canvasser groups
results.df.tj$varnum<- with(results.df.tj, paste0(as.numeric(factor(Variable))))
results.df.tj$varnum <- as.numeric(results.df.tj$varnum)
results.df.tj$sample <- rep(1:3, 8)
results.df.tj$Population <- mapvalues(results.df.tj$sample, c(1,2,3), c("All", "Left", "Right"), warn_missing = TRUE)
results.df.tj$xpos<- as.numeric(paste(results.df.tj$varnum, results.df.tj$sample, sep = "."))
results.df.tj$se.high <- results.df.tj$point.estimate + results.df.tj$se
results.df.tj$se.low <- results.df.tj$point.estimate - results.df.tj$se
results.df.tj$ci.high <- results.df.tj$point.estimate + results.df.tj$se * 1.96
results.df.tj$ci.low <- results.df.tj$point.estimate - results.df.tj$se * 1.96
labels <- unique(results.df.tj$Variable)

q <- ggplot(results.df.tj,
            aes(x=xpos, y=point.estimate,
                group=Variable, color=Population)) +
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  # CIs
  geom_linerange(aes(ymin=se.low, ymax=se.high), lwd=1) +
  geom_linerange(aes(ymin=ci.low, ymax=ci.high)) +
  # Point estimate points
  geom_point(color="black") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  ggplot2::annotate("text",
           label = c("accountable (0-4) \n ctrl mean=2.79", "advance (0-4) \n ctrl mean=1.94", "compensation (0-3) \n ctrl mean=1.93 ",
                     "force apology (0-3) \n ctrl mean=1.83", "force compensation (0-3) \n ctrl mean = 1.51", "pardoned (0-3) \n ctrl mean=0.56",
                     "public apology (0-3) \n ctrl mean=1.99", "punish (0-3) \n ctrl mean=2.19"),
           x = c(1:8)+.28, y = -1.5,
           colour = "black", size = 2.6) +
  xlab("Variable") + ylab("Coefficient") +
  ggtitle("Transitional Justice")
q +  scale_color_grey()+theme_bw()
ggsave("tj.pdf", width = 10, height = 6)

##############################################################
###### Figure 3: Positive (a) and negative (b) emotions treatment effects
##############################################################

##### Part (a) - Positive emotions
## all positive
positive <- est.ate(all$positive, all$pre_positive, all)
interested <- est.ate(all$interested, all$pre_interested, all)
stimulated <- est.ate(all$stimulated, all$pre_stimulated, all)
enthusiastic <- est.ate(all$enthusiastic, all$pre_enthusiastic, all)
energetic <- est.ate(all$energetic, all$pre_energetic, all)
proud <- est.ate(all$proud, all$pre_proud, all)
alert <- est.ate(all$alert, all$pre_alert, all)
inspired <- est.ate(all$inspired, all$pre_inspired, all)
decided <- est.ate(all$decided, all$pre_decided, all)
attentive <- est.ate(all$attentive, all$pre_attentive, all)
active <- est.ate(all$active, all$pre_active, all)

## right positive
positive.right <- est.ate(right$positive, right$pre_positive, right)
interested.right <- est.ate(right$interested, right$pre_interested, right)
stimulated.right <- est.ate(right$stimulated, right$pre_stimulated, right)
enthusiastic.right <- est.ate(right$enthusiastic, right$pre_enthusiastic, right)
energetic.right <- est.ate(right$energetic, right$pre_energetic, right)
proud.right <- est.ate(right$proud, right$pre_proud, right)
alert.right <- est.ate(right$alert, right$pre_alert, right)
inspired.right <- est.ate(right$inspired, right$pre_inspired, right)
decided.right <- est.ate(right$decided, right$pre_decided, right)
attentive.right <- est.ate(right$attentive, right$pre_attentive, right)
active.right <- est.ate(right$active, right$pre_active, right)

## left positive
positive.left <- est.ate(left$positive, left$pre_positive, left)
interested.left <- est.ate(left$interested, left$pre_interested, left)
stimulated.left <- est.ate(left$stimulated, left$pre_stimulated, left)
enthusiastic.left <- est.ate(left$enthusiastic, left$pre_enthusiastic, left)
energetic.left <- est.ate(left$energetic, left$pre_energetic, left)
proud.left <- est.ate(left$proud, left$pre_proud, left)
alert.left <- est.ate(left$alert, left$pre_alert, left)
inspired.left <- est.ate(left$inspired, left$pre_inspired, left)
decided.left <- est.ate(left$decided, left$pre_decided, left)
attentive.left <- est.ate(left$attentive, left$pre_attentive, left)
active.left <- est.ate(left$active, left$pre_active, left)

results.df.pos <- as.data.frame(rbind(positive[1], positive.left[1], positive.right[1],
                                            interested[1],interested.left[1], interested.right[1],
                                            stimulated[1], stimulated.left[1], stimulated.right[1],
                                            enthusiastic[1],  enthusiastic.left[1], enthusiastic.right[1],
                                            energetic[1], energetic.left[1],  energetic.right[1],
                                            proud[1], proud.left[1],  proud.right[1],
                                            alert[1], alert.left[1], alert.right[1],
                                            inspired[1], inspired.left[1], inspired.right[1],
                                            decided[1],  decided.left[1],decided.right[1],
                                            attentive[1], attentive.left[1],attentive.right[1], 
                                            active[1], active.left[1],active.right[1]))
results.df.pos$se <- c(positive[2], positive.left[2], positive.right[2],
                      interested[2],interested.left[2], interested.right[2],
                      stimulated[2], stimulated.left[2], stimulated.right[2],
                      enthusiastic[2],  enthusiastic.left[2], enthusiastic.right[2],
                      energetic[2], energetic.left[2],  energetic.right[2],
                      proud[2], proud.left[2],  proud.right[2],
                      alert[2], alert.left[2], alert.right[2],
                      inspired[2], inspired.left[2], inspired.right[2],
                      decided[2],  decided.left[2],decided.right[2],
                      attentive[2], attentive.left[2],attentive.right[2], 
                      active[2], active.left[2],active.right[2])

results.df.pos$se <- unlist(results.df.pos$se)
results.df.pos$point.estimate <- unlist(results.df.pos$point.estimate)
results.df.pos$Variable <- NA
results.df.pos$xpos <- NA
for (i in 1:3){results.df.pos$Variable[i] <- "positive"}
for (i in 4:6){results.df.pos$Variable[i] <- "interested"}
for (i in 7:9){results.df.pos$Variable[i] <- "stimulated"}
for (i in 10:12){results.df.pos$Variable[i] <- "enthusiastic"}
for (i in 13:15){results.df.pos$Variable[i] <- "energetic"}
for (i in 16:18){results.df.pos$Variable[i] <- "proud"}
for (i in 19:21){results.df.pos$Variable[i] <- "alert"}
for (i in 22:24){results.df.pos$Variable[i] <- "inspired"}
for (i in 25:27){results.df.pos$Variable[i] <- "decided"}
for (i in 28:30){results.df.pos$Variable[i] <- "attentive"}
for (i in 31:33){results.df.pos$Variable[i] <- "active"}

# X position of different canvasser groups
results.df.pos$varnum<- with(results.df.pos, paste0(as.numeric(factor(Variable))))
results.df.pos$varnum <- as.numeric(results.df.pos$varnum)
results.df.pos$sample <- rep(1:3, 11)
results.df.pos$Population <- mapvalues(results.df.pos$sample, c(1,2,3), c("All", "Left", "Right"), warn_missing = TRUE)
results.df.pos$xpos<- as.numeric(paste(results.df.pos$varnum, results.df.pos$sample, sep = "."))
results.df.pos$se.high <- results.df.pos$point.estimate + results.df.pos$se
results.df.pos$se.low <- results.df.pos$point.estimate - results.df.pos$se
results.df.pos$ci.high <- results.df.pos$point.estimate + results.df.pos$se * 1.96
results.df.pos$ci.low <- results.df.pos$point.estimate - results.df.pos$se * 1.96
labels <- unique(results.df.pos$Variable)
results.df.pos[1:3, 4] = results.df.pos$xpos[1:3] - 8
results.df.pos[31:33, 4] = results.df.pos$xpos[31:33] + 8

pos <- ggplot(results.df.pos,
              aes(x=xpos, y=point.estimate,
                  group=Variable, color=Population)) +
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), 
        plot.title=element_blank()) +
  # CIs
  geom_linerange(aes(ymin=se.low, ymax=se.high), lwd=1) +
  geom_linerange(aes(ymin=ci.low, ymax=ci.high)) +
  # Point estimate points
  geom_point(color="black") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  ggplot2::annotate("text",
           label = c("positive (aggregated) \n ctrl mean=18.76", "alert \n ctrl mean=1.26", "attentive \n ctrl mean=2.41", 
                     "decided \n ctrl mean=2.14", "energetic \n ctrl mean=1.81","enthusiastic \n ctrl mean=2.01",
                     "inspired \n ctrl mean=1.52", "interested \n ctrl mean=2.53",
                     "active \n ctrl mean=2.07","proud \n ctrl mean=1.24", "stimulated \n ctrl mean=1.70"),
           x = c(1:11)+.28, y = -5.8,
           colour = "black", size = 2.6) +
  xlab("Variable") + ylab("Coefficient") +ggtitle("Positive Emotions")+scale_color_grey()+theme_bw()
pos <- pos + labs(title="Positive Emotions")
ggsave("pos.pdf", width = 10, height = 6)

##### Part (b) - Negative emotions

negative <- est.ate(all$negative, all$pre_negative, all)
tense <- est.ate(all$tense, all$pre_tense, all)
scared <- est.ate(all$scared, all$pre_scared, all)
guilty <- est.ate(all$guilty, all$pre_guilty, all)
hostile <- est.ate(all$hostile, all$pre_hostile, all)
irritable <- est.ate(all$irritable, all$pre_irritable, all)
nervous <- est.ate(all$nervous, all$pre_nervous, all)
fearful <- est.ate(all$fearful, all$pre_fearful, all)
disgusted <- est.ate(all$disgusted, all$pre_disgusted, all)
afraid <- est.ate(all$afraid, all$pre_afraid, all)
embarrassed <- est.ate(all$embarrassed, all$pre_embarrassed, all)

## right negative
negative.right <- est.ate(right$negative, right$pre_negative, right)
tense.right <- est.ate(right$tense, right$pre_tense, right)
scared.right <- est.ate(right$scared, right$pre_scared, right)
guilty.right <- est.ate(right$guilty, right$pre_guilty, right)
hostile.right <- est.ate(right$hostile, right$pre_hostile, right)
irritable.right <- est.ate(right$irritable, right$pre_irritable, right)
nervous.right <- est.ate(right$nervous, right$pre_nervous, right)
fearful.right <- est.ate(right$fearful, right$pre_fearful, right)
disgusted.right <- est.ate(right$disgusted, right$pre_disgusted, right)
afraid.right <- est.ate(right$afraid, right$pre_afraid, right)
embarrassed.right <- est.ate(right$embarrassed, right$pre_embarrassed, right)

## left negative
negative.left <- est.ate(left$negative, left$pre_negative, left)
tense.left <- est.ate(left$tense, left$pre_tense, left)
scared.left <- est.ate(left$scared, left$pre_scared, left)
guilty.left <- est.ate(left$guilty, left$pre_guilty, left)
hostile.left <- est.ate(left$hostile, left$pre_hostile, left)
irritable.left <- est.ate(left$irritable, left$pre_irritable, left)
nervous.left <- est.ate(left$nervous, left$pre_nervous, left)
fearful.left <- est.ate(left$fearful, left$pre_fearful, left)
disgusted.left <- est.ate(left$disgusted, left$pre_disgusted, left)
afraid.left <- est.ate(left$afraid, left$pre_afraid, left)
embarrassed.left <- est.ate(left$embarrassed, left$pre_embarrassed, left)

# Make DF of summary stats
results.df.neg <- as.data.frame(rbind(negative[1], negative.left[1], negative.right[1],
                                            tense[1], tense.left[1], tense.right[1],
                                            scared[1],scared.left[1], scared.right[1],
                                            guilty[1], guilty.left[1], guilty.right[1],
                                            hostile[1], hostile.left[1], hostile.right[1],
                                            irritable[1], irritable.left[1], irritable.right[1],
                                            nervous[1], nervous.left[1], nervous.right[1],
                                            fearful[1],  fearful.left[1], fearful.right[1],
                                            disgusted[1], disgusted.left[1],disgusted.right[1],
                                            afraid[1], afraid.left[1],afraid.right[1],
                                            embarrassed[1], embarrassed.left[1], embarrassed.right[1]))
results.df.neg$se <- c(negative[2], negative.left[2], negative.right[2],
                                      tense[2], tense.left[2], tense.right[2],
                                      scared[2],scared.left[2], scared.right[2],
                                      guilty[2], guilty.left[2], guilty.right[2],
                                      hostile[2], hostile.left[2], hostile.right[2],
                                      irritable[2], irritable.left[2], irritable.right[2],
                                      nervous[2], nervous.left[2], nervous.right[2],
                                      fearful[2],  fearful.left[2], fearful.right[2],
                                      disgusted[2], disgusted.left[2],disgusted.right[2],
                                      afraid[2], afraid.left[2],afraid.right[2],
                                      embarrassed[2], embarrassed.left[2], embarrassed.right[2])

results.df.neg$point.estimate <- unlist(results.df.neg$point.estimate)
results.df.neg$se <- unlist(results.df.neg$se)
results.df.neg$Variable <- NA
results.df.neg$xneg <- NA
for (i in 1:3){results.df.neg$Variable[i] <- "negative"}
for (i in 4:6){results.df.neg$Variable[i] <- "tense"}
for (i in 7:9){results.df.neg$Variable[i] <- "scared"}
for (i in 10:12){results.df.neg$Variable[i] <- "guilty"}
for (i in 13:15){results.df.neg$Variable[i] <- "hostile"}
for (i in 16:18){results.df.neg$Variable[i] <- "irritable"}
for (i in 19:21){results.df.neg$Variable[i] <- "nervous"}
for (i in 22:24){results.df.neg$Variable[i] <- "fearful"}
for (i in 25:27){results.df.neg$Variable[i] <- "disgusted"}
for (i in 28:30){results.df.neg$Variable[i] <- "afraid"}
for (i in 31:33){results.df.neg$Variable[i] <- "embarrassed"}

results.df.neg$varnum<- with(results.df.neg, paste0(as.numeric(factor(Variable))))
results.df.neg$varnum <- as.numeric(results.df.neg$varnum)
results.df.neg$sample <- rep(1:3, 11)
results.df.neg$Population <- mapvalues(results.df.neg$sample, c(1,2,3), c("All", "Left", "Right"), warn_missing = TRUE)
results.df.neg$xneg<- as.numeric(paste(results.df.neg$varnum, results.df.neg$sample, sep = "."))
results.df.neg$se.high <- results.df.neg$point.estimate + results.df.neg$se
results.df.neg$se.low <- results.df.neg$point.estimate - results.df.neg$se
results.df.neg$ci.high <- results.df.neg$point.estimate + results.df.neg$se * 1.96
results.df.neg$ci.low <- results.df.neg$point.estimate - results.df.neg$se * 1.96
labels <- unique(results.df.neg$Variable)
results.df.neg[1:3, 4] = results.df.neg$xneg[1:3] - 7
results.df.neg[28:30, 4] = results.df.neg$xneg[28:30] + 7

neg <- ggplot(results.df.neg,
              aes(x=xneg, y=point.estimate,
                  group=Variable, color=Population)) +
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  # CIs
  geom_linerange(aes(ymin=se.low, ymax=se.high), lwd=1) +
  geom_linerange(aes(ymin=ci.low, ymax=ci.high)) +
  # Point estimate points
  geom_point(color="black") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  ggplot2::annotate("text",
           label = c("negative \n ctrl mean=3.44","disgusted \n ctrl mean=0.29", "embarrassed \n ctrl mean=0.26",
                     "fearful \n ctrl mean=0.22","guilty \n ctrl mean=0.24","hostile \n ctrl mean=0.33",
                     "irritable \n ctrl mean=0.54","afraid \n ctrl mean=0.16", "nervous \n ctrl mean=0.56",
                     "scared \n ctrl mean = 0.16","tense \n ctrl mean = 0.68"),
           x = c(1:11)+.28, y = -1.5,
           colour = "black", size = 2.6) +
  xlab("Variable") + ylab("Coefficient") 
neg +  scale_color_grey()+theme_bw()
ggsave("neg.pdf", width = 10, height = 6)


##############################################################
###### Figure 4: Durability of treatment effects
##############################################################

### this estimates ATE for follow up variables
est.ate.f <-function(dv, predv){
  predv <- f(predv)
  summary(fit.1 <- lm(dv~all$treat + all$pre_ideology_1 +
                        + predv*all$date_diff + all$base_gender +all$age + all$v))
  vcv <- vcovHC(fit.1)
  n <- nobs(fit.1)
  result <- coeftest(fit.1, vcv)[2, 1:4] / itt.d
  result <- list("point.estimate" = result[1],"se"=result[2],"pvalue"=result[4], "obs" = n)
  return(result)
}

# estimates ATE for follow up when we don't have a pre-treatment measurement
est.ate.np.f<-function(dv){
  summary(fit.1 <- lm(dv~all$treat + all$pre_ideology_1 + all$base_gender +all$age+all$v))
  vcv <- vcovHC(fit.1)
  n <- nobs(fit.1)
  result <- coeftest(fit.1, vcv)[2, 1:4] / itt.d
  result <- list("point.estimate" = result[1],"se"=result[2],"pvalue"=result[4], "obs" = n)
  return(result)
}

par(mfrow=c(4, 2))
par(oma = c(2.9, 3, 1, 0)) # make room (i.e. the 4's) for the overall x and y axis titles
par(mar = c(2.5, 2, 1.5, 1)) # make the plots be closer together
par(cex.main = 0.8)

##### GRAPH 1 - CHURCH TRUST ######

church_trust_pre <- est.ate.np(all$pre_conf_church, all)
church_trust <- est.ate(all$conf_church, all$pre_conf_church,all)
church_trust_f1 <- est.ate.f(all$conf_church_f1, all$pre_conf_church)
church_trust_f2 <- est.ate.f(all$conf_church_f2, all$pre_conf_church)
church_trust_f3 <- est.ate.f(all$conf_church_f3, all$pre_conf_church)

coefs <- unlist(c(church_trust_pre[1], church_trust[1], church_trust_f1[1],
           church_trust_f2[1], church_trust_f3[1]))

ses <- unlist(c(church_trust_pre[2], church_trust[2], church_trust_f1[2],
         church_trust_f2[2], church_trust_f3[2]))


plot(NA, xlim = c(-2, 25), ylim = c(-.5, .5),xlab = '', ylab = '')
title("Trust in church (0-3)")
abline(v = -1, col = "gray")
abline(v = 0, col = "gray")
abline(v = 1, col = "gray")
abline(v = 8, col = "gray")
abline(v = 24, col = "gray")
abline(h = 0, col = "red")

points(-1, coefs[1], pch = 23, col = "black", bg = "black")
points(0, coefs[2], pch = 23, col = "black", bg = "black")
points(1, coefs[3], pch = 23, col = "black", bg = "black")
points(8, coefs[4], pch = 23, col = "black", bg = "black")
points(24, coefs[5], pch = 23, col = "black", bg = "black")

segments(-1, (coefs - ses)[1], -1, (coefs + ses)[1], col = "black", lwd = 2)
segments(0, (coefs - ses)[2], 0, (coefs + ses)[2], col = "black", lwd = 2)
segments(1, (coefs - ses)[3], 1, (coefs + ses)[3], col = "black", lwd = 2)
segments(8.0, (coefs - ses)[4], 8, (coefs + ses)[4], col = "black", lwd = 2)
segments(24.0, (coefs - ses)[5], 24, (coefs + ses)[5], col = "black", lwd = 2)
segments(-1, (coefs - 1.96*ses)[1], -1, (coefs + 1.96*ses)[1], col = "black", lwd = 1)
segments(0, (coefs - 1.96*ses)[2], 0, (coefs + 1.96*ses)[2], col = "black", lwd = 1)
segments(1, (coefs - 1.96*ses)[3], 1, (coefs + 1.96*ses)[3], col = "black", lwd = 1)
segments(8.0, (coefs - 1.96*ses)[4], 8, (coefs + 1.96*ses)[4], col = "black", lwd = 1)
segments(24.0, (coefs - 1.96*ses)[5], 24, (coefs + 1.96*ses)[5], col = "black", lwd = 1)

text(-2, -.2, "pre", cex = .6)
text(-1, -.25, "treatment", cex = .6)
text(.2, .3, "treatment", cex = .6)
text(3, .07, "follow up", cex = .6)
text(10, .07, "follow up", cex = .6)
text(22, .13, "follow up", cex = .6)

##### GRAPH 2 - PARDONING ######

pardoned <- est.ate.np(all$policies_pardon,all)
pardoned_f1<- est.ate.np.f(all$policies_pardon_f1)
pardoned_f2 <- est.ate.np.f(all$policies_pardon_f2)
pardoned_f3 <- est.ate.np.f(all$policies_pardon_f3)

coefs <- unlist(c(pardoned[1], pardoned_f1[1],
           pardoned_f2[1], pardoned_f3[1]))

ses <- unlist(c(pardoned[2], pardoned_f1[2],
         pardoned_f2[2], pardoned_f3[2]))

plot(NA, xlim = c(-.5, 25), ylim = c(-.2, .5), xlab = '', ylab = '')
title("Support for pardoning perpetrators (0-4)")
abline(v = 0, col = "gray")
abline(v = 1, col = "gray")
abline(v = 8, col = "gray")
abline(v = 24, col = "gray")
abline(h = 0, col = "red")

points(0, coefs[1], pch = 23, col = "black", bg = "black")
points(1, coefs[2], pch = 23, col = "black", bg = "black")
points(8, coefs[3], pch = 23, col = "black", bg = "black")
points(24, coefs[4], pch = 23, col = "black", bg = "black")

segments(0, (coefs - ses)[1], 0, (coefs + ses)[1], col = "black", lwd = 2)
segments(1, (coefs - ses)[2], 1, (coefs + ses)[2], col = "black", lwd = 2)
segments(8.0, (coefs - ses)[3], 8, (coefs + ses)[3], col = "black", lwd = 2)
segments(24.0, (coefs - ses)[4], 24, (coefs + ses)[4], col = "black", lwd = 2)

segments(0, (coefs - 1.96*ses)[1], 0, (coefs + 1.96*ses)[1], col = "black", lwd = 1)
segments(1, (coefs - 1.96*ses)[2], 1, (coefs + 1.96*ses)[2], col = "black", lwd = 1)
segments(8.0, (coefs - 1.96*ses)[3], 8, (coefs + 1.96*ses)[3], col = "black", lwd = 1)
segments(24.0, (coefs - 1.96*ses)[4], 24, (coefs + 1.96*ses)[4], col = "black", lwd = 1)

##### GRAPH 3 - GOVT SATISFACTION ######
gov_sat_pre <- est.ate.np(all$pre_inst_gov,all)
gov_sat <- est.ate(all$inst_gov, all$pre_inst_gov,all)
gov_sat_f1 <- est.ate.f(all$inst_gov_f1, all$pre_inst_gov)
gov_sat_f2 <- est.ate.f(all$inst_gov_f2, all$pre_inst_gov)
gov_sat_f3 <- est.ate.f(all$inst_gov_f3, all$pre_inst_gov)

coefs <- unlist(c(gov_sat_pre[1], gov_sat[1], gov_sat_f1[1],
           gov_sat_f2[1], gov_sat_f3[1]))

ses <- unlist(c(gov_sat_pre[2], gov_sat[2], gov_sat_f1[2],
         gov_sat_f2[2], gov_sat_f3[2]))


plot(NA, xlim = c(-2, 25), ylim = c(-.4, .5), xlab = '', ylab = '')
title("Satisfaction with government (0-3)")
abline(v = -1, col = "gray")
abline(v = 0, col = "gray")
abline(v = 1, col = "gray")
abline(v = 8, col = "gray")
abline(v = 24, col = "gray")
abline(h = 0, col = "red")

points(-1, coefs[1], pch = 23, col = "black", bg = "black")
points(0, coefs[2], pch = 23, col = "black", bg = "black")
points(1, coefs[3], pch = 23, col = "black", bg = "black")
points(8, coefs[4], pch = 23, col = "black", bg = "black")
points(24, coefs[5], pch = 23, col = "black", bg = "black")

segments(-1, (coefs - ses)[1], -1, (coefs + ses)[1], col = "black", lwd = 2)
segments(0, (coefs - ses)[2], 0, (coefs + ses)[2], col = "black", lwd = 2)
segments(1, (coefs - ses)[3], 1, (coefs + ses)[3], col = "black", lwd = 2)
segments(8.0, (coefs - ses)[4], 8, (coefs + ses)[4], col = "black", lwd = 2)
segments(24.0, (coefs - ses)[5], 24, (coefs + ses)[5], col = "black", lwd = 2)
segments(-1, (coefs - 1.96*ses)[1], -1, (coefs + 1.96*ses)[1], col = "black", lwd = 1)
segments(0, (coefs - 1.96*ses)[2], 0, (coefs + 1.96*ses)[2], col = "black", lwd = 1)
segments(1, (coefs - 1.96*ses)[3], 1, (coefs + 1.96*ses)[3], col = "black", lwd = 1)
segments(8.0, (coefs - 1.96*ses)[4], 8, (coefs + 1.96*ses)[4], col = "black", lwd = 1)
segments(24.0, (coefs - 1.96*ses)[5], 24, (coefs + 1.96*ses)[5], col = "black", lwd = 1)

##### GRAPH 4 - DEMOCRACY ######
dem_pre <- est.ate.np(all$pre_democracy,all)
dem <- est.ate(all$democracy, all$pre_democracy,all)
dem_f1 <- est.ate.f(all$democracy_f1, all$pre_democracy)
dem_f2 <- est.ate.f(all$democracy_f2, all$pre_democracy)
dem_f3 <- est.ate.f(all$democracy_f3, all$pre_democracy)

coefs <- unlist(c(dem_pre[1], dem[1], dem_f1[1],dem_f2[1], dem_f3[1]))

ses <- unlist(c(dem_pre[2], dem[2], dem_f1[2],dem_f2[2], dem_f3[2]))

plot(NA, xlim = c(-2, 25), ylim = c(-.4, .3), xlab = '', ylab = '')
title("Satisfaction with democracy (0-3)")
abline(v = -1, col = "gray")
abline(v = 0, col = "gray")
abline(v = 1, col = "gray")
abline(v = 8, col = "gray")
abline(v = 24, col = "gray")
abline(h = 0, col = "red")

points(-1, coefs[1], pch = 23, col = "black", bg = "black")
points(0, coefs[2], pch = 23, col = "black", bg = "black")
points(1, coefs[3], pch = 23, col = "black", bg = "black")
points(8, coefs[4], pch = 23, col = "black", bg = "black")
points(24, coefs[5], pch = 23, col = "black", bg = "black")

segments(-1, (coefs - ses)[1], -1, (coefs + ses)[1], col = "black", lwd = 2)
segments(0, (coefs - ses)[2], 0, (coefs + ses)[2], col = "black", lwd = 2)
segments(1, (coefs - ses)[3], 1, (coefs + ses)[3], col = "black", lwd = 2)
segments(8.0, (coefs - ses)[4], 8, (coefs + ses)[4], col = "black", lwd = 2)
segments(24.0, (coefs - ses)[5], 24, (coefs + ses)[5], col = "black", lwd = 2)
segments(-1, (coefs - 1.96*ses)[1], -1, (coefs + 1.96*ses)[1], col = "black", lwd = 1)
segments(0, (coefs - 1.96*ses)[2], 0, (coefs + 1.96*ses)[2], col = "black", lwd = 1)
segments(1, (coefs - 1.96*ses)[3], 1, (coefs + 1.96*ses)[3], col = "black", lwd = 1)
segments(8.0, (coefs - 1.96*ses)[4], 8, (coefs + 1.96*ses)[4], col = "black", lwd = 1)
segments(24.0, (coefs - 1.96*ses)[5], 24, (coefs + 1.96*ses)[5], col = "black", lwd = 1)

##### GRAPH 5 - MIL GOV ######
mil_pre <- est.ate.np(all$pre_military_gov,all)
mil <- est.ate(all$military_gov, all$pre_military_gov,all)
mil_f1 <- est.ate.f(all$military_gov_f1, all$pre_military_gov)
mil_f2 <- est.ate.f(all$military_gov_f2, all$pre_military_gov)
mil_f3 <- est.ate.f(all$military_gov_f3, all$pre_military_gov)

coefs <- unlist(c(mil_pre[1], mil[1], mil_f1[1],mil_f2[1], mil_f3[1]))
ses <- unlist(c(mil_pre[2], mil[2], mil_f1[2],mil_f2[2], mil_f3[2]))

plot(NA, xlim = c(-2, 25), ylim = c(-.2, .3), xlab = '', ylab = '')
title("Support for military gov (0-1)")
abline(v = -1, col = "gray")
abline(v = 0, col = "gray")
abline(v = 1, col = "gray")
abline(v = 8, col = "gray")
abline(v = 24, col = "gray")
abline(h = 0, col = "red")

points(-1, coefs[1], pch = 23, col = "black", bg = "black")
points(0, coefs[2], pch = 23, col = "black", bg = "black")
points(1, coefs[3], pch = 23, col = "black", bg = "black")
points(8, coefs[4], pch = 23, col = "black", bg = "black")
points(24, coefs[5], pch = 23, col = "black", bg = "black")

segments(-1, (coefs - ses)[1], -1, (coefs + ses)[1], col = "black", lwd = 2)
segments(0, (coefs - ses)[2], 0, (coefs + ses)[2], col = "black", lwd = 2)
segments(1, (coefs - ses)[3], 1, (coefs + ses)[3], col = "black", lwd = 2)
segments(8.0, (coefs - ses)[4], 8, (coefs + ses)[4], col = "black", lwd = 2)
segments(24.0, (coefs - ses)[5], 24, (coefs + ses)[5], col = "black", lwd = 2)

segments(-1, (coefs - 1.96*ses)[1], -1, (coefs + 1.96*ses)[1], col = "black", lwd = 1)
segments(0, (coefs - 1.96*ses)[2], 0, (coefs + 1.96*ses)[2], col = "black", lwd = 1)
segments(1, (coefs - 1.96*ses)[3], 1, (coefs + 1.96*ses)[3], col = "black", lwd = 1)
segments(8.0, (coefs - 1.96*ses)[4], 8, (coefs + 1.96*ses)[4], col = "black", lwd = 1)
segments(24.0, (coefs - 1.96*ses)[5], 24, (coefs + 1.96*ses)[5], col = "black", lwd = 1)

##### GRAPH 6 - ADVANCE ######
advance <- est.ate.np(all$justice_advance,all)
advance_f1 <- est.ate.np.f(all$justice_advance_f1)
advance_f2 <- est.ate.np.f(all$justice_advance_f2)
advance_f3 <- est.ate.np.f(all$justice_advance_f3)

coefs <- unlist(c(advance[1], advance_f1[1], advance_f2[1], advance_f3[1]))

ses <- unlist(c(advance[2], advance_f1[2], advance_f2[2], advance_f3[2]))

plot(NA, xlim = c(-.5, 25), ylim = c(-.8, .2), xlab = '', ylab = '')
title("Obsession with the past \n makes it hard to advance (0-4)")
abline(v = 0, col = "gray")
abline(v = 1, col = "gray")
abline(v = 8, col = "gray")
abline(v = 24, col = "gray")
abline(h = 0, col = "red")

points(0, coefs[1], pch = 23, col = "black", bg = "black")
points(1, coefs[2], pch = 23, col = "black", bg = "black")
points(8, coefs[3], pch = 23, col = "black", bg = "black")
points(24, coefs[4], pch = 23, col = "black", bg = "black")

segments(0, (coefs - ses)[1], 0, (coefs + ses)[1], col = "black", lwd = 2)
segments(1, (coefs - ses)[2], 1, (coefs + ses)[2], col = "black", lwd = 2)
segments(8.0, (coefs - ses)[3], 8, (coefs + ses)[3], col = "black", lwd = 2)
segments(24, (coefs - ses)[4], 24, (coefs + ses)[4], col = "black", lwd = 2)
segments(0, (coefs - 1.96*ses)[1], 0, (coefs + 1.96*ses)[1], col = "black", lwd = 1)
segments(1, (coefs - 1.96*ses)[2], 1, (coefs + 1.96*ses)[2], col = "black", lwd = 1)
segments(8.0, (coefs - 1.96*ses)[3], 8, (coefs + 1.96*ses)[3], col = "black", lwd = 1)
segments(24.0, (coefs - 1.96*ses)[4], 24, (coefs + 1.96*ses)[4], col = "black", lwd = 1)

## GRAPH 7 - COMPENSATION ##
comp_pre <- est.ate.np(all$pre_current_recomp,all)
comp <- est.ate(all$current_recomp, all$pre_current_recomp,all)
comp_f1 <- est.ate.f(all$current_recomp_f1, all$pre_current_recomp)
comp_f2 <- est.ate.f(all$current_recomp_f2, all$pre_current_recomp)
comp_f3 <- est.ate.f(all$current_recomp_f3, all$pre_current_recomp)

coefs <- unlist(c(comp_pre[1], comp[1], comp_f1[1],comp_f2[1], comp_f3[1]))

ses <- unlist(c(comp_pre[2], comp[2], comp_f1[2], comp_f2[2], comp_f3[2]))

plot(NA, xlim = c(-2, 25), ylim = c(-.4, .4), xlab = '', ylab = '')
title("Support for victim compensation (0-3)")
abline(v = -1, col = "gray")
abline(v = 0, col = "gray")
abline(v = 1, col = "gray")
abline(v = 8, col = "gray")
abline(v = 24, col = "gray")
abline(h = 0, col = "red")

points(-1, coefs[1], pch = 23, col = "black", bg = "black")
points(0, coefs[2], pch = 23, col = "black", bg = "black")
points(1, coefs[3], pch = 23, col = "black", bg = "black")
points(8, coefs[4], pch = 23, col = "black", bg = "black")
points(24, coefs[5], pch = 23, col = "black", bg = "black")

segments(-1, (coefs - ses)[1], -1, (coefs + ses)[1], col = "black", lwd = 2)
segments(0, (coefs - ses)[2], 0, (coefs + ses)[2], col = "black", lwd = 1)
segments(1, (coefs - ses)[3], 1, (coefs + ses)[3], col = "black", lwd = 2)
segments(8.0, (coefs - ses)[4], 8, (coefs + ses)[4], col = "black", lwd = 2)
segments(24.0, (coefs - ses)[5], 24, (coefs + ses)[5], col = "black", lwd = 2)
segments(-1, (coefs - 1.96*ses)[1], -1, (coefs + 1.96*ses)[1], col = "black", lwd = 1)
segments(0, (coefs - 1.96*ses)[2], 0, (coefs + 1.96*ses)[2], col = "black", lwd = 1)
segments(1, (coefs - 1.96*ses)[3], 1, (coefs + 1.96*ses)[3], col = "black", lwd = 1)
segments(8.0, (coefs - 1.96*ses)[4], 8, (coefs + 1.96*ses)[4], col = "black", lwd = 1)
segments(24.0, (coefs - 1.96*ses)[5], 24, (coefs + 1.96*ses)[5], col = "black", lwd = 1)

## GRAPH 8 - PUBLIC APOLOGY ##

tj_apology <- est.ate.np(all$tj_apology,all)
tj_apology_f1 <- est.ate.np.f(all$tj_apology_f1)
tj_apology_f2 <- est.ate.np.f(all$tj_apology_f2)
tj_apology_f3 <- est.ate.np.f(all$tj_apology_f3)

coefs <- unlist(c(tj_apology[1], tj_apology_f1[1], tj_apology_f2[1], tj_apology_f3[1]))

ses <- unlist(c(tj_apology[2], tj_apology_f1[2], tj_apology_f2[2], tj_apology_f3[2]))

plot(NA, xlim = c(-.5, 25), ylim = c(-.2, .5), xlab = '', ylab = '')
title("Military should apologize (0-4)")
abline(v = 0, col = "gray")
abline(v = 1, col = "gray")
abline(v = 8, col = "gray")
abline(v = 24, col = "gray")
abline(h = 0, col = "red")

points(0, coefs[1], pch = 23, col = "black", bg = "black")
points(1, coefs[2], pch = 23, col = "black", bg = "black")
points(8, coefs[3], pch = 23, col = "black", bg = "black")
points(24, coefs[4], pch = 23, col = "black", bg = "black")

segments(0, (coefs - ses)[1], 0, (coefs + ses)[1], col = "black", lwd = 2)
segments(1, (coefs - ses)[2], 1, (coefs + ses)[2], col = "black", lwd = 2)
segments(8.0, (coefs - ses)[3], 8, (coefs + ses)[3], col = "black", lwd = 2)

segments(24.0, (coefs - ses)[4], 24, (coefs + ses)[4], col = "black", lwd = 2)
segments(0, (coefs - 1.96*ses)[1], 0, (coefs + 1.96*ses)[1], col = "black", lwd = 1)
segments(1, (coefs - 1.96*ses)[2], 1, (coefs + 1.96*ses)[2], col = "black", lwd = 1)
segments(8.0, (coefs - 1.96*ses)[3], 8, (coefs + 1.96*ses)[3], col = "black", lwd = 1)
segments(24.0, (coefs - 1.96*ses)[4], 24, (coefs + 1.96*ses)[4], col = "black", lwd = 1)

mtext('Weeks from Treatment', side = 1, outer = TRUE, line = 2)
mtext('Coefficient', side = 2, outer = TRUE, line = 2)